home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / swagg_m.zip / MATH.SWG / 0016_PI2.PAS.pas < prev    next >
Pascal/Delphi Source File  |  1993-05-28  |  3KB  |  190 lines
  1. {
  2.         Here's a good (but a little slow) Program to calculate the
  3.   decimals of Pi :
  4.  
  5.  
  6. THIS Program CompUTES THE DIGITS of PI USinG THE ARCTANGENT ForMULA
  7. (1)            PI/4 = 4 ARCTAN 1/5 - ARCTAN 1/239
  8. in CONJUNCTION With THE GREGorY SERIES
  9.  
  10. (2)   ARCTAN X = SUM  (-1)^N*(2N + 1)^-1*X^(2N+1)  N=0 to  inFinITY.
  11.  
  12. SUBSTITUTinG into (2) A FEW VALUES of N  and NESTinG  WE HAVE,
  13.  
  14. PI/4 =  1/5[4/1 + 1/25[-4/3 + 1/25[4/5 + 1/25[-4/7 + ...].].]
  15.  
  16.     - 1/239[1/1 + 1/239^2[-1/3 + 1/239^2[1/5 + 1/239^2[-1/7 +...].].]
  17.  
  18. USinG THE LONG divISION ALGorITHM, THIS ( NESTED ) inFinITE SERIES CAN BE
  19. USED to CALCULATE PI to A LARGE NUMBER of DECIMAL PLACES in A REASONABLE
  20. AMOUNT of TIME. A TIME Function IS inCLUDED to SHOW HOW SLOW THinGS
  21. GET WHEN N IS LARGE. IMPROVEMENTS CAN BE MADE BY CHANGinG THE SIZE of
  22. THE Array ELEMENTS HOWEVER IT GETS A BIT TRICKY.
  23.  
  24. }
  25.  
  26. Uses
  27.   Crt;
  28.  
  29. Var
  30.   B,C,V,P1,S,K,N,I,J,Q,M,M1,X,R,D : Integer;
  31.   P,A,T : Array[0..5000] of Integer;
  32.  
  33. Const F1=5;
  34. Const F2=239;
  35. Procedure divIDE(D : Integer);
  36.  begin
  37.     R:=0;
  38.     For J:=0 to M do
  39.      begin
  40.      V:= R*10+P[J];
  41.      Q:=V div D;
  42.      R:=V Mod D;
  43.      P[J]:=Q;
  44.      end;
  45. end;
  46. Procedure divIDEA(D : Integer);
  47.  begin
  48.     R:=0;
  49.     For J:=0 to M do
  50.      begin
  51.      V:= R*10+A[J];
  52.      Q:=V div D;
  53.      R:=V Mod D;
  54.      A[J]:=Q;
  55.      end;
  56.  end;
  57. Procedure SUBT;
  58. begin
  59. B:=0;
  60. For J:=M Downto 0 do
  61.     if T[J]>=A[J]  then T[J]:=T[J]-A[J] else
  62.     begin
  63.      T[J]:=10+T[J]-A[J];
  64.      T[J-1]:=T[J-1]-1;
  65.    end;
  66. For J:=0 to M do
  67. A[J]:=T[J];
  68. end;
  69. Procedure SUBA;
  70. begin
  71. For J:=M Downto 0 do
  72.     if P[J]>=A[J]  then P[J]:=P[J]-A[J] else
  73.     begin
  74.      P[J]:=10+P[J]-A[J];
  75.      P[J-1]:=P[J-1]-1;
  76.    end;
  77. For J:= M Downto 0 do
  78. A[J]:=P[J];
  79. end;
  80. Procedure CLEARP;
  81.  begin
  82.   For J:=0 to M do
  83.    P[J]:=0;
  84.  end;
  85. Procedure ADJUST;
  86. begin
  87. P[0]:=3;
  88. P[M]:=10;
  89. For J:=1 to M-1 do
  90. P[J]:=9;
  91. end;
  92.  
  93. Procedure ADJUST2;
  94. begin
  95. P[0]:=0;
  96. P[M]:=10;
  97. For J:=1 to M-1 do
  98. P[J]:=9;
  99. end;
  100.  
  101. Procedure MULT4;
  102.  begin
  103.   C:=0;
  104.   For J:=M Downto 0 do
  105.    begin
  106.     P1:=4*A[J]+C;
  107.     A[J]:=P1 Mod 10;
  108.     C:=P1 div 10;
  109.    end;
  110.   end;
  111.  
  112. Procedure SAVEA;
  113. begin
  114. For J:=0 to M do
  115. T[J]:=A[J];
  116. end;
  117.  
  118. Procedure TERM1;
  119. begin
  120.  I:=M+M+1;
  121.  A[0]:=4;
  122. divIDEA(I*25);
  123. While I>3 do
  124. begin
  125. I:=I-2;
  126. CLEARP;
  127. P[0]:=4;
  128. divIDE(I);
  129. SUBA;
  130. divIDEA(25);
  131. end;
  132. CLEARP;
  133. ADJUST;
  134. SUBA;
  135. divIDEA(5);
  136. SAVEA;
  137. end;
  138. Procedure TERM2;
  139. begin
  140.  I:=M+M+1;
  141.  A[0]:=1;
  142. divIDEA(I);
  143. divIDEA(239);
  144. divIDEA(239);
  145. While I>3 do
  146. begin
  147. I:=I-2;
  148. CLEARP;
  149. P[0]:=1;
  150. divIDE(I);
  151. SUBA;
  152. divIDEA(239);
  153. divIDEA(239);
  154. end;
  155. CLEARP;
  156. ADJUST2;
  157. SUBA;
  158. divIDEA(239);
  159. SUBT;
  160. end;
  161.  
  162. {MAin Program}
  163.  
  164.    begin
  165.    ClrScr;
  166.    WriteLn('                        THE CompUTATION of PI');
  167.    WriteLn('                     -----------------------------');
  168.    WriteLn('inPUT NO. DECIMAL PLACES');
  169.    READLN(M1);
  170.    M:=M1+4;
  171.     For J:=0 to M  do
  172.        begin
  173.          A[J]:=0;
  174.          T[J]:=0;
  175.        end;
  176.    TERM1;
  177.    TERM2;
  178.    MULT4;
  179.    WriteLn;WriteLn;
  180.    Write('PI = 3.');
  181.    For J:=1 to M1   do
  182.    begin
  183.     Write(A[J]);
  184.    if J Mod 5 =0 then Write(' ');
  185.    if J Mod 50=0 then Write('                    ');
  186.    end;
  187.    WriteLn;WriteLn;
  188.    WriteLn;
  189. end.
  190.